// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

#include "Tetrahedralizer.h"
#include "mat3d_tet.h"
#include <vtkCell.h>
#include <VisItException.h>

#include <algorithm>
#include <vector>
using std::vector;

// ****************************************************************************
//  Constructor: Tetrahedralizer::Tetrahedralizer
//
//  Programmer:  Jeremy Meredith
//  Creation:    August 30, 2002
//
// ****************************************************************************
Tetrahedralizer::Tetrahedralizer(int nmat_)
{
    nmat = nmat_;
    for (int i=0; i<MAX_TETS_PER_CELL; i++)
        tet[i] = MaterialTetrahedron(nmat);
}

// ****************************************************************************
//  Method:  Tetrahedralizer::GetLowTetNodesForTet
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a tetrahedron.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    December 11, 2000
//
//  Modifications:
//    Jeremy Meredith, Tue Sep 18 12:01:19 PDT 2001
//    Added 'origindex' so we can avoid using the weights when possible.
//
//    Jeremy Meredith, Tue Oct 16 09:05:20 PDT 2001
//    Added a clear() before the resize so the mattet's get re-initialized
//    every time.  This lets us re-use the Tetrahedralizer across cells.
//
//    Jeremy Meredith, Fri Aug 30 17:03:45 PDT 2002
//    Removed resizing of vector, and instead reset the output tets.
//
//    Jeremy Meredith, Tue Jan 14 15:03:22 PST 2003
//    Extracted to a public static class function and made it fill in an
//    array instead of creating all the tets, so that it could be used
//    in more contexts.
//
// ****************************************************************************
int
Tetrahedralizer::GetLowTetNodesForTet(int nnodes,const vtkIdType *nodes, vtkIdType *tetnodes)
{
    int ntet = 1;
    tetnodes[0] = 0;
    tetnodes[1] = 1;
    tetnodes[2] = 2;
    tetnodes[3] = 3;
    return ntet;
}

// ****************************************************************************
//  Method:  Tetrahedralizer::GetLowTetNodesForVox
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a voxel.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    December 11, 2000
//
//  Modifications:
//    Jeremy Meredith, Tue Sep 18 12:01:19 PDT 2001
//    Added 'origindex' so we can avoid using the weights when possible.
//
//    Jeremy Meredith, Tue Oct 16 09:05:20 PDT 2001
//    Added a clear() before the resize so the mattet's get re-initialized
//    every time.  This lets us re-use the Tetrahedralizer across cells.
//
//    Jeremy Meredith, Fri Aug 30 17:03:45 PDT 2002
//    Removed resizing of vector, and instead reset the output tets.
//
//    Jeremy Meredith, Mon Jan  6 10:07:21 PST 2003
//    VTK tets are inverted from Silo tets, so I flipped the ordering of
//    all created tets.
//
//    Jeremy Meredith, Tue Jan 14 15:03:22 PST 2003
//    Extracted to a public static class function and made it fill in an
//    array instead of creating all the tets, so that it could be used
//    in more contexts.
//
// ****************************************************************************
int
Tetrahedralizer::GetLowTetNodesForVox(int nnodes,const vtkIdType *nodes, vtkIdType *tetnodes)
{
    int n;

    // find the node with the lowest index
    int minv = nodes[0];
    int minn = 0;
    for (n=1; n<nnodes; n++)
    {
        if (minv > nodes[n])
        {
            minv = nodes[n];
            minn = n;
        }
    }

    // rearrange the zone so that node is now at the zero position
    int flipx = (minn==0 || minn==4 || minn==2 || minn==6) ? 0 : 1;
    int flipy = (minn==4 || minn==5 || minn==6 || minn==7) ? 0 : 1;
    int flipz = (minn==0 || minn==1 || minn==4 || minn==5) ? 0 : 1;
    
    // create the old-node to new-node cross-index
    int *map  = mat3d_vox2tet_subdiv1_flip[flipx][flipy][flipz];

    // after flipping, is the zone going to be inside-out?
    bool invert = ((flipx + flipy + flipz) % 2) == 0;

    // which diagonals will the tets cross...
    int diagx0 = std::min(nodes[map[1]],nodes[map[6]]);
    int diagx1 = std::min(nodes[map[2]],nodes[map[5]]);
    int diagy0 = std::min(nodes[map[1]],nodes[map[3]]);
    int diagy1 = std::min(nodes[map[2]],nodes[map[0]]);
    int diagz0 = std::min(nodes[map[6]],nodes[map[3]]);
    int diagz1 = std::min(nodes[map[5]],nodes[map[0]]);

    int diagx = (diagx0 < diagx1) ? 0 : 1;
    int diagy = (diagy0 < diagy1) ? 0 : 1;
    int diagz = (diagz0 < diagz1) ? 0 : 1;

    // ... and which case do these correspond to (after flipping)
    int tetcase = (diagx<<2) + (diagy<<1) + (diagz);
    
    int ntet = ntet_mat3d_vox2tet_subdiv1[tetcase];

    for (int t=0; t<ntet; t++)
    {
        if (!invert)
        {
            tetnodes[t*4+0] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][2]];
            tetnodes[t*4+3] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][3]];
        }
        else
        {
            tetnodes[t*4+0] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][3]];
            tetnodes[t*4+3] = map[mat3d_vox2tet_subdiv1_node[tetcase][t][2]];
        }
    }
    return ntet;
}

// ****************************************************************************
//  Method:  Tetrahedralizer::GetLowTetNodesForHex
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a hexahedron.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    December 11, 2000
//
//  Modifications:
//    Jeremy Meredith, Mon Feb 12 17:11:36 PST 2001
//    Fixed a numbering problem causing adjacent faces to be tetrahedralized
//    differently.
//
//    Jeremy Meredith, Tue Sep 18 12:01:19 PDT 2001
//    Added 'origindex' so we can avoid using the weights when possible.
//
//    Jeremy Meredith, Tue Oct 16 09:05:20 PDT 2001
//    Added a clear() before the resize so the mattet's get re-initialized
//    every time.  This lets us re-use the Tetrahedralizer across cells.
//
//    Jeremy Meredith, Fri Aug 30 17:03:45 PDT 2002
//    Removed resizing of vector, and instead reset the output tets.
//
//    Jeremy Meredith, Mon Jan  6 10:07:21 PST 2003
//    VTK tets are inverted from Silo tets, so I flipped the ordering of
//    all created tets.
//
//    Jeremy Meredith, Tue Jan 14 15:03:22 PST 2003
//    Extracted to a public static class function and made it fill in an
//    array instead of creating all the tets, so that it could be used
//    in more contexts.
//
// ****************************************************************************
int
Tetrahedralizer::GetLowTetNodesForHex(int nnodes,const vtkIdType *nodes, vtkIdType *tetnodes)
{
    int n;

    // find the node with the lowest index
    int minv = nodes[0];
    int minn = 0;
    for (n=1; n<nnodes; n++)
    {
        if (minv > nodes[n])
        {
            minv = nodes[n];
            minn = n;
        }
    }

    // rearrange the zone so that node is now at the zero position
    int flipx = (minn==0 || minn==4 || minn==7 || minn==3) ? 0 : 1;
    int flipy = (minn==5 || minn==4 || minn==7 || minn==6) ? 0 : 1;
    int flipz = (minn==0 || minn==1 || minn==4 || minn==5) ? 0 : 1;
    
    // create the old-node to new-node cross-index
    int *map  = mat3d_hex2tet_subdiv1_flip[flipx][flipy][flipz];

    // after flipping, is the zone going to be inside-out?
    bool invert = ((flipx + flipy + flipz) % 2) == 0;

    // which diagonals will the tets cross...
    int diagx0 = std::min(nodes[map[1]],nodes[map[6]]);
    int diagx1 = std::min(nodes[map[2]],nodes[map[5]]);
    int diagy0 = std::min(nodes[map[1]],nodes[map[3]]);
    int diagy1 = std::min(nodes[map[2]],nodes[map[0]]);
    int diagz0 = std::min(nodes[map[6]],nodes[map[3]]);
    int diagz1 = std::min(nodes[map[2]],nodes[map[7]]);

    int diagx = (diagx0 < diagx1) ? 0 : 1;
    int diagy = (diagy0 < diagy1) ? 0 : 1;
    int diagz = (diagz0 < diagz1) ? 0 : 1;

    // ... and which case do these correspond to (after flipping)
    int tetcase = (diagx<<2) + (diagy<<1) + (diagz);
    
    int ntet = ntet_mat3d_hex2tet_subdiv1[tetcase];

    for (int t=0; t<ntet; t++)
    {
        if (!invert)
        {
            tetnodes[t*4+0] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][2]];
            tetnodes[t*4+3] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][3]];
        }
        else
        {
            tetnodes[t*4+0] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][3]];
            tetnodes[t*4+3] = map[mat3d_hex2tet_subdiv1_node[tetcase][t][2]];
        }
    }
    return ntet;
}

// ****************************************************************************
//  Method:  Tetrahedralizer::GetLowTetNodesForWdg
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a wedge.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    December 11, 2000
//
//  Modifications:
//    Jeremy Meredith, Tue Sep 18 12:01:19 PDT 2001
//    Added 'origindex' so we can avoid using the weights when possible.
//
//    Jeremy Meredith, Tue Oct 16 09:05:20 PDT 2001
//    Added a clear() before the resize so the mattet's get re-initialized
//    every time.  This lets us re-use the Tetrahedralizer across cells.
//
//    Jeremy Meredith, Fri Aug 30 17:03:45 PDT 2002
//    Removed resizing of vector, and instead reset the output tets.
//
//    Jeremy Meredith, Mon Jan  6 10:07:21 PST 2003
//    VTK tets are inverted from Silo tets, so I flipped the ordering of
//    all created tets.
//
//    Jeremy Meredith, Tue Jan 14 15:03:22 PST 2003
//    Extracted to a public static class function and made it fill in an
//    array instead of creating all the tets, so that it could be used
//    in more contexts.
//
// ****************************************************************************
int
Tetrahedralizer::GetLowTetNodesForWdg(int nnodes,const vtkIdType *nodes, vtkIdType *tetnodes)
{
    int n;

    // find the node with the lowest index
    int minv = nodes[0];
    int minn = 0;
    for (n=1; n<nnodes; n++)
    {
        if (minv > nodes[n])
        {
            minv = nodes[n];
            minn = n;
        }
    }

    // rearrange the zone so that node is now at the zero position
    int flip  = (minn==0 || minn==1 || minn==2) ? 0 : 1;
    int rot   = (minn==0 || minn==3) ? 0 : ((minn==1 || minn==4) ? 1 : 2);
    
    // create the old-node to new-node cross-index
    int *map  = mat3d_wedge2tet_subdiv1_fliprot[flip][rot];

    // after flipping, is the zone going to be inside-out?
    bool invert = (flip % 2) == 0;

    // which diagonals will the tets cross...
    int diag0  = std::min(nodes[map[2]],nodes[map[4]]);
    int diag1  = std::min(nodes[map[1]],nodes[map[5]]);

    int diag   = (diag0 < diag1) ? 0 : 1;

    // ... and which case does this correspond to (after flipping)
    int tetcase = diag;
    
    int ntet = ntet_mat3d_wedge2tet_subdiv1;

    for (int t=0; t<ntet; t++)
    {
        if (!invert)
        {
            tetnodes[t*4+0] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][2]];
            tetnodes[t*4+3] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][3]];
        }
        else
        {
            tetnodes[t*4+0] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][3]];
            tetnodes[t*4+3] = map[mat3d_wedge2tet_subdiv1_node[tetcase][t][2]];
        }
    }
    return ntet;
}

// ****************************************************************************
//  Method:  Tetrahedralizer::GetLowTetNodesForPyr
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a pyramid.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    December 11, 2000
//
//  Modifications:
//    Jeremy Meredith, Tue Sep 18 12:01:19 PDT 2001
//    Added 'origindex' so we can avoid using the weights when possible.
//
//    Jeremy Meredith, Tue Oct 16 09:05:20 PDT 2001
//    Added a clear() before the resize so the mattet's get re-initialized
//    every time.  This lets us re-use the Tetrahedralizer across cells.
//
//    Jeremy Meredith, Fri Aug 30 17:03:45 PDT 2002
//    Removed resizing of vector, and instead reset the output tets.
//
//    Jeremy Meredith, Mon Jan  6 10:07:21 PST 2003
//    VTK tets are inverted from Silo tets, so I flipped the ordering of
//    all created tets.
//
//    Jeremy Meredith, Tue Jan 14 15:03:22 PST 2003
//    Extracted to a public static class function and made it fill in an
//    array instead of creating all the tets, so that it could be used
//    in more contexts.
//
// ****************************************************************************
int
Tetrahedralizer::GetLowTetNodesForPyr(int nnodes,const vtkIdType *nodes, vtkIdType *tetnodes)
{
    int n;

    // find the node with the lowest index
    vtkIdType minv = nodes[0];
    int minn = 0;
    for (n=1; n<nnodes; n++)
    {
        if (minv > nodes[n])
        {
            minv = nodes[n];
            minn = n;
        }
    }

    // rearrange the zone so that node is now at the zero position
    int flip  = (minn==0 || minn==2) ? 0 : 1;
    
    // create the old-node to new-node cross-index
    int *map  = mat3d_pyr2tet_subdiv1_flip[flip];

    // after flipping, is the zone going to be inside-out?
    bool invert = (flip % 2) == 0;

    // which diagonals will the tets cross...
    // -- no question for a pyramid
    // ... and which case this must correspond (after flipping)
    int tetcase = 0;
    
    int ntet = ntet_mat3d_pyr2tet_subdiv1;

    for (int t=0; t<ntet; t++)
    {
        if (!invert)
        {
            tetnodes[t*4+0] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][2]];
            tetnodes[t*4+3] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][3]];
        }
        else
        {
            tetnodes[t*4+0] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][0]];
            tetnodes[t*4+1] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][1]];
            tetnodes[t*4+2] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][3]];
            tetnodes[t*4+3] = map[mat3d_pyr2tet_subdiv1_node[tetcase][t][2]];
        }
    }
    return ntet;
}


// ****************************************************************************
//  Method:  Tetrahedralizer::calc_tet_low
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a tetrahedron and
//    fill in the tets.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    January 14, 2003
//
//  Modifications:
//
// ****************************************************************************
void
Tetrahedralizer::calc_tet_low()
{
    ntet=GetLowTetNodesForTet(nnodes,nodes,tetnodes_low);
    fill_tets_low();
}

// ****************************************************************************
//  Method:  Tetrahedralizer::calc_vox_low
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a voxel and
//    fill in the tets.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    January 14, 2003
//
//  Modifications:
//
// ****************************************************************************
void
Tetrahedralizer::calc_vox_low()
{
    ntet=GetLowTetNodesForVox(nnodes,nodes,tetnodes_low);
    fill_tets_low();
}

// ****************************************************************************
//  Method:  Tetrahedralizer::calc_hex_low
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a hexahedron and
//    fill in the tets.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    January 14, 2003
//
//  Modifications:
//
// ****************************************************************************
void
Tetrahedralizer::calc_hex_low()
{
    ntet=GetLowTetNodesForHex(nnodes,nodes,tetnodes_low);
    fill_tets_low();
}

// ****************************************************************************
//  Method:  Tetrahedralizer::calc_wdg_low
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a wedge and
//    fill in the tets.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    January 14, 2003
//
//  Modifications:
//
// ****************************************************************************
void
Tetrahedralizer::calc_wdg_low()
{
    ntet=GetLowTetNodesForWdg(nnodes,nodes,tetnodes_low);
    fill_tets_low();
}

// ****************************************************************************
//  Method:  Tetrahedralizer::calc_pyr_low
//
//  Purpose:
//    Calculate the low-subdivision tetrahedralization for a pyramid and
//    fill in the tets.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    January 14, 2003
//
//  Modifications:
//
// ****************************************************************************
void
Tetrahedralizer::calc_pyr_low()
{
    ntet=GetLowTetNodesForPyr(nnodes,nodes,tetnodes_low);    
    fill_tets_low();
}

// ****************************************************************************
//  Method:  fill_tets_low
//
//  Purpose:
//    Fills in the tetrahedron array from the precalculated
//    ntet & tetnodes_low arrays.
//
//  Arguments:
//    
//
//  Programmer:  Jeremy Meredith
//  Creation:    January 14, 2003
//
// ****************************************************************************
void
Tetrahedralizer::fill_tets_low()
{
    for (int t=0; t<ntet; t++)
    {
        tet[t].reset();

        int a = tetnodes_low[t*4+0];
        int b = tetnodes_low[t*4+1];
        int c = tetnodes_low[t*4+2];
        int d = tetnodes_low[t*4+3];

        tet[t].node[0].origindex = nodes[a];
        tet[t].node[1].origindex = nodes[b];
        tet[t].node[2].origindex = nodes[c];
        tet[t].node[3].origindex = nodes[d];

        tet[t].node[0].weight[a] = 1.0;
        tet[t].node[1].weight[b] = 1.0;
        tet[t].node[2].weight[c] = 1.0;
        tet[t].node[3].weight[d] = 1.0;

        for (int m=0; m<nmat; m++)
        {
            tet[t].node[0].matvf[m] = vf_node[a][m];
            tet[t].node[1].matvf[m] = vf_node[b][m];
            tet[t].node[2].matvf[m] = vf_node[c][m];
            tet[t].node[3].matvf[m] = vf_node[d][m];
        }
    }
}

// ****************************************************************************
//  Method:  Tetrahedralizer::Tetrahedralize
//
//  Purpose:
//    Create the tetrahedralization using the given parameters.
//
//  Arguments:
//    s          the subivision level
//    c          the cell
//    nmat_      the number of materials
//    vf_xxxx_   the locally indexed VFs for the zone/node/face/edge
//
//  Programmer:  Jeremy Meredith
//  Creation:    December 11, 2000
//
//  Modifications:
//    Hank Childs, Mon Oct  7 17:43:41 PDT 2002 
//    Remove costly VTK calls.
//
//    Hank Childs, Fri Jan 28 15:43:04 PST 2005
//    Use exception macros.
//
// ****************************************************************************
void
Tetrahedralizer::Tetrahedralize(MIROptions::SubdivisionLevel s, int celltype,
                                int npts, const vtkIdType *c_ptr,
                                vector<float>   vf_zone_,
                                vector<float>  *vf_node_,
                                vector<float>  *vf_face_,
                                vector<float>  *vf_edge_)
{
    subdiv = s;
    nnodes = npts;
    nodes  = c_ptr;

    vf_zone = vf_zone_;
    vf_node = vf_node_;
    vf_face = vf_face_;
    vf_edge = vf_edge_;

    switch (subdiv)
    {
      case MIROptions::Low:
        switch (celltype)
        {
          case VTK_TETRA:       calc_tet_low();   break;
          case VTK_VOXEL:       calc_vox_low();   break;
          case VTK_HEXAHEDRON:  calc_hex_low();   break;
          case VTK_WEDGE:       calc_wdg_low();   break;
          case VTK_PYRAMID:     calc_pyr_low();   break;
          default:              EXCEPTION0(VisItException);
        }
        break;
      case MIROptions::Med:
        //
        // If you implement subdivision level medium you must change
        // the code in MIR::ReconstructTet.
        //
        switch (celltype)
        {
          case VTK_TETRA:       ntet = 12;     break;
          case VTK_VOXEL:       ntet = 24;     break;
          case VTK_HEXAHEDRON:  ntet = 24;     break;
          case VTK_WEDGE:       ntet = 18;     break;
          case VTK_PYRAMID:     ntet = 13;     break;
          default:              EXCEPTION0(VisItException);
        }
        break;
      case MIROptions::High:
        //
        // If you implement subdivision level high you must change
        // the code in MIR::ReconstructTet.
        //
        switch (celltype)
        {
          case VTK_TETRA:       ntet = 24;     break;
          case VTK_VOXEL:       ntet = 48;     break;
          case VTK_HEXAHEDRON:  ntet = 48;     break;
          case VTK_WEDGE:       ntet = 36;     break;
          case VTK_PYRAMID:     ntet = 26;     break;
          default:              EXCEPTION0(VisItException);
        }
        break;
      default:
        EXCEPTION0(VisItException);
    }
}
